Exploiting Large-Scale Correlations to Detect Continuous Gravitational Waves 
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Fully coherent searches (over realistic ranges of parameter space and year-long observation times) for 
unknown sources of continuous gravitational waves are computationally prohibitive. Less expensive hi- 
erarchical searches divide the data into shorter segments which are analyzed coherently, then detection 
statistics from different segments are combined incoherently. The novel method presented here solves the 
long-standing problem of how best to do the incoherent combination. The optimal solution exploits large- 
scale parameter-space correlations in the coherent detection statistic. Application to simulated data shows 
dramatic sensitivity improvements compared with previously available (ad hoc) methods, increasing the 
spatial volume probed by more than 2 orders of magnitude at lower computational cost. 
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Searching for CW Sources. — Direct detection of 
gravitational waves is the most significant remaining test 
of Einstein's General Theory of Relativity, and will become 
an important new astronomical tool. 

Rapidly rotating neutron stars are expected to gener- 
ate continuous gravitational- wave (CW) signals via various 
mechanisms 11] |2l [3 IH [51 . Most such stars are electro- 
magnetically invisible, but might be detected and studied 
via gravitational waves. Recent simulations of neutron star 
populations [6, 7| suggest that CW sources might eventu- 
ally be detected with new instruments such as LIGO [8, 9|. 
World-wide efforts are underway to search for CW sig- 
nals llOl nil [12j and observational upper limits already 
place some constraints on neutron star physics [13, 14J. 

Because the expected CW signals are weak, sensitive 
data analysis methods are needed to extract these sig- 
nals from detector noise. A powerful method is derived 
in Ref. [15|. This scheme is based on the principle of 
maximum likelihood detection, which leads to coherent 
matched filtering. Rotating neutron stars emit monochro- 
matic CW signals, apart from a slowly changing intrinsic 
frequency. But the terrestrial detector location Doppler- 
modulates the amplitude and phase of the waveform, as the 
Earth moves relative to the solar system barycenter (SSB). 
The parameters describing the signal's amplitude variation 
may be analytically eliminated by maximizing the coherent 
matched-filtering statistic [15J. The remaining search pa- 
rameters describing the signal's phase are the source's sky 
location, frequency and frequency derivatives. The result- 
ing coherent detection statistic is called the JT-statistic. 

This work considers isolated CW sources whose fre- 
quency varies linearly with time in the SSB frame. 
The corresponding phase parameter- space V is four- 
dimensional. Standard "physical" coordinates on V are 
the frequency f(to) at some fiducial time to, the fre- 
quency's first time derivative /, and a unit vector n = 
(cos 5 cos a, cos 5 sin a, sin 5) on the two-sphere 5^, 
pointing from the SSB to the source. Here a and S are 
right ascension and declination. Thus, a point in parameter 



space p G P may be labeled by p = {f{to), /, n}. The 
JT-statistic [h] is a functional of the detector data set h, 
and is a function of the point in parameter space p G 7^. 

All-sky searches for unknown CW sources using the 
jT-statistic are computationally expensive. For maximum 
sensitivity, one must convolve the full data set with sig- 
nal waveforms (templates) corresponding to all possible 
sources. But the number of templates required for a fully 
coherent search increases as a high power of the observa- 
tion time. For one year of data, the computational cost to 
search a realistic range of parameter space exceeds the total 
computing power on Earth |[T5][T6l. Thus a fully coherent 
search is limited to much shorter observation times. 

Searching year-long data sets is accomplished by less 
costly hierarchical semicoherent methods |[T7l[T8l[T9l. The 
data is broken into segments of duration T, where T is 
much smaller than one year. Each segment is analyzed 
coherently, computing the JT-statistic on a coarse grid of 
templates. Then the values from all segments (or statis- 
tics derived from JT) are incoherently combined using a 
common fine grid of templates, discarding phase informa- 
tion between segments. Among the current semicoherent 
strategies, the Stack-Slide method [17, 18 ] sums values 
along putative signal tracks in the time-frequency plane. 
The Hough transform method fT9l sums H{!F — T^h) 
where JT^h is a constant predefined threshold. The Heavy- 
side function i?(x) is unity for positive x and vanishes 
elsewhere. This latter technique is currently used by Ein- 
stein© Home [12J, a public distributed computing project 
carrying out the most sensitive all-sky CW searches. 

A central long-standing problem in these semicoherent 
methods is the design of, and link between, the coarse and 
fine grids. Current methods, while creative and clever, are 
arbitrary and ad hoc constructions. This work removes 
all arbitrariness by finding the optimal solution for the 
incoherent step. The key quantity is the fractional loss, 
called mismatch Ai, in expected JT-statistic (or sum of 
jT-statistics in the incoherent step) for a given signal p 
at a nearby grid point p^ Locally Taylor-expanding A4 
(to quadratic order) in the differences of the coordinates 



{/(to), / , n} of p and defines a signature ++++ metric 
ds^ [[T6l|20l[2Tll221- Current methods consider parameter 
correlations in T to only linear order in T . 

The JT-statistic has large-scale correlations ||23l [231 in 
the physical coordinates {/(to), /, n}, extending outside 
a region in which the mismatch is well-approximated by 
the metric given above ll25l . Recent work ll24l has shown 
that (for a given signal) the region where the expected T- 
statistic has maximal value may be described by a separate 
equation for each order of T, when T is small compared 
to one year. The solutions to these equations are hyper- 
surfaces in V , whose intersections are the extrema of (an 
approximation to) T . 

For currently used values of T (a day or longer) it is also 
crucial to consider the fractional loss of T to second order 
in T [24J. For source frequencies > 1 kHz and for values 
of T > 60 h, additional orders in T would be needed. 

The new method. — This work exploits the large-scale 
correlations in the coherent detection statistic T to con- 
struct a significantly improved semicoherent search tech- 
nique for CW signals. The novel method is optimal if the 
semicoherent detection statistic is taken to be the sum of 
one coarse-grid JT-statistic value from each data segment, 
and makes four important improvements. 

First, the improved understanding of large-scale correla- 
tions yields new coordinates on V. The metric in these co- 
ordinates accurately approximates the mismatch in T [ 25 1 
in each segment. Hence, the optimal (closest) coarse-grid 
point from each segment can be determined for any given 
fine-grid point in the incoherent combination step. 

Second, in the new coordinates we find the first analyt- 
ical metric for the semicoherent detection statistic to con- 
struct the optimal fine grid (minimum possible number of 
points). Previous ad hoc approaches obtain the ne grid by 
rening the coarse grid in three dimensions, / and n. Here, 
the explicit incoherent- step metric shows that refinement is 
only needed in one dimension, /. This greatly reduces the 
computational cost at equal detection sensitivity, although 
it also reduces the accuracy of the estimated source param- 
eters. But this is a very profitable trade, because in a hi- 
erarchical search the primary goal of the first stages is to 
discard the uninteresting regions of parameter space. Later 
follow-up stages use longer coherent integrations to more 
accurately determine the source parameters. 

Third, existing techniques combine the coherent results 
less effectively than our method, because they do not use 
metric information beyond linear order in T. This leads to 
a higher sensitivity at equal computational cost. 

Fourth, the new technique can simultaneously do a 
Stack-Slide-like summing of values and a Hough-like 
summing of i7(jF — JF^h) , with a lower total computational 
cost than either one of these methods individually. 

For a given CW source with realistic phase parameter 
values (/ ^ 1 kHz, |/| < //50 yr) and coherent data seg- 
ment lengths T < 60 h, the large-scale correlations of the 



jT-statistic are well described by the first- and second-order 
(in T) equations |24|: 

K^) = /(^) + /(^)i(^)-n + /|(t)-n, 

^(^) = / + /(^)l(^)-n + 2/i(t)-n, 

where/(t) = /(to) + (t-to)/. (1) 

Here ^(t) = rorb(^)/c, with rorb(^) denoting the vector 
from the Earth's barycenter to the SSB, and c the speed of 
light. The quantities jy{t) and z>(t) can be interpreted as the 
source's instantaneous frequency and frequency derivative 
at the Earth's barycenter at time t. 

The parameters u and z> provide new coordinates on V. 
It is also useful to introduce new (real- valued) sky coordi- 
nates rix and riy (as in [26]): 

nx(t) +zny(t) = f{t)rE cos 6d cos^e^t^-^^^'^l (2) 

Here Te = Re/c ^ 21ms is the light travel time from 
the Earth's center to the detector, and a^it), are the 
detector position at time t. The metric separation ds'^ is 

d^V^' =dz/' TVs + 7' dz>' rVl80 + 2 dn^ + 2 dn^ 
- 4 dz/ duy T/{7t£) + 4 dz> drix T^/{i^lf . 

(3) 

In defining differences in coordinates {z/, z>, rix, riy}, the 
time t in Eqs. ([T]) and ([2]) is the midpoint of the data seg- 
ment spanning times [t — T/2, t + T/2], and 7 = 1. To 
simplify the form of the metric, T is taken to be a positive 
integer number I of sidereal days. 

The new coordinates {z/, z>, rix, riy} have important ad- 
vantages over the original coordinates {/, /, n}. The met- 
ric is explicitly coordinate-independent (showing that V is 
flat). In fact, the region around a point p in which the mis- 
match M is well-approximated by ds^ is much larger [251 . 

Consider a segment of data which contains a strong 
CW source with phase parameters p. If the sky separation 
patch is small enough to neglect the drix and dUy terms in 
Eq. ([3]), then jFp/ [h^] is extremized for all that have the 
same values of v and z> as p. This set of points in V forms 
a two-dimensional surface dv — dv — 0. Thus, for all 
sources within the sky patch, there exists a different (/, /) 
pair with those same values of v and z>. This property is 
exploited by the new method. 

An implementation. — The data set is divided into 
N segments of length T (potentially including short gaps 
in operation) labeled by the integer j = 1, A^. The 
segments span time intervals [tj — T/2, + r/2]. The 
detector-time midpoint of segment j is tj and to = 

X^jLi is fiducial time. 

Every segment is analyzed coherently on a coarse grid 
in V. This grid is constructed so that no point in V is 
farther than a specified distance from some coarse-grid 
point, where the distance measure is defined by the met- 
ric of Eq. ([3]). To simplify the grid construction, large fre- 
quency bands are analyzed by breaking them into many 
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narrow sub-bands. For each segment j, and at each coarse 
grid point, the JF-statistic is evaluated, and "stats" are 
obtained. Here, the word "stat" denotes the two-tuple 
{Tj,H{Tj - J^th)). 

For simplicity, the same coarse grid is used for all data 
segments: the Cartesian product of a rectangular grid in 
/, / and a grid on the sky-sphere n G 5^. The spac- 
ings are A/ = Vvim/{TiT) and A/ = V7mm/{7vT^), 
where m is the one-dimensional metric mismatch param- 
eter inn . The spacing of the coarse sky grid is chosen so 
that the dn^ and dUy terms in Eq. ^ may be neglected. 
When orthogonally projected onto the equatorial unit disk, 
the sky grid is uniform, and contains ^ 27r/(A(/^)^ points, 
with A(f = V2m/ (tt/ Te cos 6^)' 

The incoherent step combines the "stats" obtained by the 
coherent analysis, using a fine grid in P. At each point in 
the fine grid, a "stat" value is obtained by summing one 
Stat value from each of the N coarse grids. The coarse 
grid point is the one with the same sky position as the fine 
grid point, which has the smallest metric separation from 
Eq. ([3]). The final result is a "stat" value at each point on 
the fine grid. The first element of the stat is the sum of the 
jT-statistic values from the coarse-grid points. The second 
element is a number count, reflecting the number of data 
segments in which jFth was exceeded. A detectable CW 
signal leads to a fine-grid point with a high number count 
and a large sum of JT-statistics. 

The spacing of the fine grid is determined from the met- 
ric for the fractional loss of the expected X^jLi 
parameter offsets between a putative signal location and a 
fine grid point. This may be calculated as proposed in [ 17 1, 
by averaging the coarse-grid metric over the N different 
segments. Since each coarse-grid metric is no longer cal- 
culated at the data-segment midpoints (but at to), the co- 
efficients change between segments because of the time- 
dependence of the parameter-space coordinates. For our 
choices of to and T, the only additional term in the metric 
Eq. ^ that does not average to zero is {tj — to)^r^dz>^/3. 
The averaged metric is the diagonal part of Eq. ([3]) with 
72 _ 1 + mYj^^^itj - toy/{NT^), where the coordi- 
nate offsets are calculated at the fiducial time to. Thus, the 
fine grid may be identical to the coarse grid except that the 
spacing A/ is smaller by a factor 7. This is of order N 
when the number of data segments is large. No further re- 
finement in frequency or sky position is needed. Coherent 
integration over the total observation time would require 
refining both A/ and A/ (increasing points oc N^), plus 
similar sky refinements. 

Performance comparison. — Monte Carlo simula- 
tions are used to illustrate the improved performance of 
this method compared with the conventional Hough trans- 
form technique [19|. The software tools used are part of 
LALApps ll^ 1 and employ accurate barycentering routines 
with timing errors below 3/xs. To provide a realistic com- 
parison, simulated data covered the same time intervals as 



the input data used for the current (S5R5) Einstein@Home 
search [12J. Those data, from LIGO Hanford (HI, 4km) 
and LIGO Livingston (LI, 4km), are not contiguous, but 
contain gaps when the detectors are not operating. The to- 
tal time interval spanned is about 264 days, containing 121 
data segments of duration 25 h (so approximately £ = 1). 

The false alarm probabilities are obtained using 5 000 
simulated data sets with different realizations of stationary 
Gaussian white noise, with one-sided strain spectral den- 
sity ^/S^ = 3.25 X 10-22 Hz"^/2 To find the detection 
probabilities, different CW signals with fixed strain am- 
plitude ho are added. The parameters [15 ] are randomly 
drawn from uniform distributions in cos(inclination l), po- 
larization i/j, initial phase 00, the entire sky, /(to) G 
[100.1, 100.3] Hz, and / G [-1.29, -0.711] nHz/s. 

Figure [T] compares the performance of the different 
methods. The receiver operating characteristic is the de- 
tection probability as a function of false alarm probability, 
at fixed source strain amplitude /lo = 6 x lO"^^^ Because 
the number count (using J^^h — 2.6) is discrete, the two 
"curves" in Fig. [T] consist of discrete points. Our method 
(using either number counts or summed as a detection 
statistic) is superior to the conventional Hough technique. 

In addition, this method is computationally faster. The 
comparison used identical coherent stages (m = 0.3, with 
2 981 coarse-grid points) for both this method and the con- 
ventional Hough. But using different fine grids in the inco- 
herent step, this method's fine grid had 506 times as many 
points as the coarse grid, but the Hough fine grid had 7 056 
times as many points. In spite of using 14 times fewer fine- 
grid points, our method has substantially higher sensitivity. 

Figure [2] compares the methods. It shows the detec- 
tion efficiencies for different values of source strain am- 
plitude /lo, at a fixed 1% false alarm probability. As above, 
each point in Fig. [2] is obtained by analyzing 2 000 sim- 
ulated data sets. Again, this technique in both modes 
of operation performs substantially better than the Hough 
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FIG. 1: Receiver operating characteristic curve. The new method 
performs better than the conventional Hough technique. 



3 




Source Strain Amplitude / 10 ' 



FIG. 2: Detection probability as a function of source strain am- 
plitude /lo, at 1% false alarm probability. 

method. For example, the ho value needed to obtain 90% 
detection probability with our technique in number-count 
operation mode is smaller by a factor ^ 6 than the ho 
needed by the Hough method: the "distance reach" ifTSll of 
our technique is ^ 6 times larger. This increases the num- 
ber of potentially detectable sources by more than 2 orders 
of magnitude, since the "visible" spatial volume increases 
as the cube of the distance. The lower computational cost 
of our method would also allow increases in N or T, even 
further improving the sensitivity. 

These results are qualitatively independent of frequency, 
as confirmed in additional comparisons. 

Conclusions. — A novel semicoherent technique for 
detecting CW sources has been described. In contrast 
to previous approaches, the new method exploits large- 
scale parameter-space correlations in the coherent detec- 
tion statistic to optimally solve the subsequent inco- 
herent combination step. For coherent integration times 
T < 60 h, the correlations are well-described by the 
second-order (in T) formulae presented here. The method 
should be extendible to longer coherent integration times 
by including higher orders. It could be extended to search 
for CW sources in binary systems, as well as to space- 
based detectors. The method also has applicability in radio, 
X- and 7-ray astronomy (searches for weak radio or 7-ray 
pulsars, or pulsations from X-ray binaries). 

Realistic Monte Carlo simulations show that our tech- 
nique is much more sensitive than the conventional Hough 
method (current most sensitive all-sky CW search tech- 
nique). The technique presented here is also computation- 
ally simpler, and more efficient. 

The LIGO Scientific Collaboration is currently work- 
ing to deploy this technique on the Einstein @ Home 
project [12J, starting with LIGO S6 data. The combina- 
tion of a better search technique, and more sensitive data, 
greatly increases the chance of making the first gravita- 



tional wave detection of a CW source. The detection of 
CW signals will provide new means to discover and locate 
neutron stars, and will eventually provide unique insights 
into the nature of matter at high densities. 

We thank R. Prix, M. A. Papa, C. Messenger, B. Kr- 
ishnan and B. Knispel for helpful discussions. This work 
was supported by the Max Planck Gesellschaft and by 
U.S. National Science Foundation Grant No. 0555655 and 
No. 0701817. 



* Electronic address: Holger.Pletsch@aei.mpg.de| 
^ Electronic address: Bruce.Allen@aei.mpg.de 
[1] B. J. Owen et al, Phys. Rev. D 58, 084020 (1998). 
[2] G. Ushomirsky, C. Cutler, and L. Bildsten, Mon. Not. Roy. 

Astron. Soc. 319, 902 (2000). 
[3] C. Cutler, Phys. Rev. D 66, 084025 (2002). 
[4] D. 1. Jones and N. Andersson, Mon. Not. Roy. Astron. Soc. 

331, 203 (2002). 
[5] B. J. Owen, Phys. Rev. Lett. 95, 211101 (2005). 
[6] B. Knispel and B. Allen, Phys. Rev. D 78, 044031 (2008). 
[7] C. J. Horowitz and K. Kadau, Phys. Rev. Lett. 102, 191102 

(2009). 

[8] A. Abramovici et al. Science 256, 325 (1992). 
[9] B. C. Barish and R. Weiss, Physics Today 52, 44 (1999). 
[10] B. Abbott et al. (The LIGO Scientific Collaboration), 

Phys. Rev D 77, 022001 (2008). 
[11] B. Abbott et al. (The LIGO Scientific Collaboration), 

Phys. Rev D 79, 022001 (2009). 
[12] Einstein @ Home: http://einstein.phys.uwm.edu/. 
[13] B. Abbott et al. (The LIGO Scientific Collaboration), Astro- 

phys. J. Lett. 683, L45 (2008). 
[14] B. Abbott et al. (The LIGO Scientific Collaboration), 

Phys. Rev Lett. 102, 111102 (2009). 
[15] P Jaranowski, A. Krolak, and B. F. Schutz, Phys. Rev D 58, 

063001 (1998). 
[16] P R. Brady et al, Phys. Rev D 57, 2101 (1998). 
[17] P R. Brady and T. Creighton, Phys. Rev D 61, 082001 

(2000). 

[18] C Cutler, 1. Gholami, and B. Krishnan, Phys. Rev. D 72, 
042004 (2005). 

[19] B. Krishnan et al., Phys. Rev D 70, 082001 (2004). 

[20] R. Balasubramanian, B. S. Sathyaprakash, and S. V. Dhu- 
randhar, Phys. Rev D 53, 3033 (1996). 

[21] B. J. Owen, Phys. Rev D 53, 6749 (1996). 

[22] R. Prix, Phys. Rev D 75, 023004 (2007). 

[23] R. Prix and Y. Itoh, Class. Quant. Grav 22, S1003 (2005). 

[24] H. J. Pletsch, Phys. Rev D 78, 102005 (2008). 

[25] Using {zy, z>, rix, riy} coordinates the metric ds^ remains ac- 
curate up to = 0.3. In contrast, in the {/, /, n} co- 
ordinates ds^ can yield errors greater than 10% for mis- 
matches as small as 0.001. In part, this is because the metric 
in {/, /, n} varies significantly over the M = 0.3 region. 

[26] P Astone et al., Phys. Rev D 65, 042003 (2002). 

[27] http://www.lsc-group.phys.uwm.edu/daswg/. 



4 



